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Abstract — Many cell functions are accomplished thanks to 
intracellular transport mechanisms of macromolecules along 
filaments. Molecular motors such as dynein or kinesin are 
proteins playing a primary role in these processes. The behavior 
of such proteins is quite well understood when there is only one 
of them moving a cargo particle. Indeed, numerous in vitro 
experiments have been performed to derive accurate models 
for a single molecular motor. However, in vivo macromolecules 
are often carried by multiple motors. The main focus of 
this paper is to provide an analysis of the behavior of more 
molecular motors interacting together in order to improve the 
understanding of their actual physiological behavior. Previous 
studies provide analyses based on results obtained from Monte 
Carlo simulations. Different from these studies, we derive 
an equipollent probabilistic model to describe the dynamics 
of multiple proteins coupled together and provide an exact 
theoretical analysis. We are capable of obtaining the probability 
density function of the motor protein configurations, thus 
enabling a deeper understanding of their behavior. 

I. Introduction 

Many biological cell functions are accomplished thanks 
to intracellular transport processes. A key role in these 
processes is played by molecular motors that are proteins 
converting chemical energy into mechanical work. Among 
the most important molecular motors, there are kinesin and 
dynein. These are classes of molecular motors transporting 
many kind of different molecules by moving along filaments 
called microtubuli. They accomplish this by converting the 
chemical energy stored in the ATP molecules produce me- 
chanical work and thus represent powerful nanomachines. 
Both kinesin and dynein have a similar structure consisting 
of two catalytically active motor domains, called "heads" 
that can bind both to a microtubule and ATP. The two heads 
are hinged together by a polypeptidic linkage that terminates 
in a distal cargo binding domain. These motor heads bind 
and unbind to the microtubule successively by utilizing the 
energy available from ATP hydrolysis and enable the motor 
to walk [1]. 

The study of molecular motors has developed over the 
past decade and these investigations have unraveled many 
previously unknown finer details at the microscopic level 
which have propelled various theoretical discussions of the 
mechanisms underlying the dynamics of molecular motors. 
Given the small sizes of molecular motors (each head of 
kinesin has a linear size of the order of 10 nm), they are 
supposed to be perturbed by the state of the local molecular 



environment and thermal fluctuations around them. But the 
directed walk of these cytoskeletal motors clearly reveals 
that they are robust enough to overcome the fluctuations 
surrounding them. The behavior of motor proteins is quite 
well understood when there is only one of them moving 
a cargo particle. Numerous in vitro experiments have been 
performed to derive accurate models to describe the motion 
of a single molecular motor. Indeed, it is possible to monitor 
the motion of a single molecular motor under highly tunable 
experimental conditions and obtain measurements with accu- 
rate spatial and time resolution [2], [3], [4]. As a result, there 
are many theoretical descriptions of motor-protein mecha- 
nisms that recognize their multiple conformational transitions 
taking into account the complex mechanochemical processes 
involved [5]. 

However, the chemical and mechanical properties of molec- 
ular motors in vivo are not yet completely characterized. 
Furthermore, there is evidence that, in vivo, macromolecules 
are carried by multiple motors [6]. Indeed, this can be 
a cause of the robustness of their behavior and of also 
a cause of increased efficiency. The possible dynamics of 
an ensemble of motor proteins carrying the same cargo 
particle is still largely unknown. For example, it is not 
clear whether synchronization phenomena occur, the motors 
move independently or a "tug-of-war" is in action [7], [8]. 
Understanding these mechanisms is one of the more chal- 
lenging problems that require concerted efforts by chemists, 
physicists, biologists and engineers. 

II. Markov model for molecular motors 

In this section we derive a Markov model to describe 
the motion of an ensemble of molecular motors carrying a 
cargo. While the behavior of single molecular motors such as 
kinesin or dynein is quite well understood, we intend to focus 
on the possible properties of synchronization/coordination 
that can be exhibited by more motors connected to the same 
cargo. 

Molecular motors are proteins moving macro-molecules and 
organelles along filaments (i.e. microtubuli). Their motion 
occurs by discrete steps. Their heads move forward by 
hydrolyzing ATP and binding to proper locations that are 
equally spaced on the microtubule. Every motor can attach 
to the cargo molecule via a flexible linkage (see Figure [TJ. 




Fig. 1. A molecular motor making a step 



In [6] a stochastic algorithm simulating the dynamics of 
multiple interacting motors was employed to obtain useful 
insights about their behavior via Monte Carlo simulations. 
We "borrow" the rules of the algorithm given by [6] and 
introduce a rigorous mathematical formalism in order to 
derive a Markov model and investigate how multiple motors 
function together. However, in contrast to [6], our approach 
will have the advantage of relying on a probabilistic analysis 
that is exact and is not based on Monte Carlo simulations. 
Every motor is modeled as a particle moving on a line by 
discrete steps. Its position can assume values that are equally 
spaced by the unit step size d s . Multiple motors are allowed 
to bind to the same location. Every motor linkage behaves 
as a spring with elastic constant K and rest length 1$ with 
no compressional rigidity: when compressed it buckles down 
with no resistance. It is assumed that the motors are bound 
to the same cargo particle. The cargo is subject to the forces 
applied by the motors and to a constant load Fi oa d that is 
assumed positive when it opposes the motor movement. 
In [6], two different models are provided for the cargo mo- 
tion. In a simplified model the cargo moves deterministically 
reaching its equilibrium position immediately. In a more 
realistic model the cargo is subject to a Brownian force that 
models the thermal fluctuations of the environment. 
In this paper we focus on the first scenario only, even though 
our results can be extended to treat the noisy case, as well. 
The motors move by single steps of length d s in a stochastic 
manner. Backward steps are not allowed. Given two non- 
overlapping time intervals, the probabilities that a motor 
takes a step in each of them are independent. This leads 
to a Poisson distribution of the time after which every next 
step occurs. For each motor, the stepping rate of the Poisson 
distribution depends on the reaction force F it is subject to 
from the cargo. F is assumed positive when it opposes the 
motor movement. A motor subject to a force F will have a 
stepping rate given by \( step '(F). Estimates for \( ste P)(F) 
can be obtained experimentally by single motor experiments 
and from the knowledge of the chemical process behind the 
ATP hydrolyzation. 

Given the forces acting on each motor, they step indepen- 
dently from each other. The motors have a "stalling force" 
F s such that, for F > F s , we have \ ste p(F) = 0. The 
stalling force can be determined via experiments, as well. 
The force applied to a motor also influence the probability 
rate X det (F) at which it can detach form the microtubule. 
Every unengaged motor can also reattach to the microtubule 
at a rate X att , but only at locations that are distant not more 
than Iq from the cargo position. In other words, a motor can 



not attach to the microtubule at a location where its linkage 
is stretched. Finally, the model assumes that at any time there 
can not be more than m motors engaged on the microtubule. 

A. Mathematical formalism and Notation 

We model the microtubule locations as a bi-infinite se- 
quence {ak}kez> a k € R> where each a k represents the 
position of a binding site at which a motor can be located. 
Since the locations are equally spaced by the distance d s , 
we provide the following definition 

Definition 1: We define a location sequence as A := 
{a k }k&, where 



ak 



r (o) 



kd s 



(1) 



and x^ E K. 

Note that x^ € K is determined by the arbitrary choice of 
the reference point on the microtubule. 
The arrangement of the motors on the microtubule is de- 
scribed by a bi-infinite sequence of natural numbers {zk}k£Z, 
Zk E N, where z k represents the number of motors bound at 
the position a&. Equivalently, the arrangement of motors on 
the microtubule can be described by a finite ordered sequence 
{xi)YL\ representing the position of each of the m motors 
engaged on the microtubule. 

Definition 2: Assume a location sequence A = {ak]k^z- 
We refer to a configuration as any bi-infinite sequence Z := 
{zk}kez of natural numbers such that 



E 



z k < oo. 



(2) 



The natural number m is the number of engaged agents of 
the configuration Z. 

We define the position sequence of the configuration Z 
with respect to the location sequence A as the unique non- 
decreasing finite sequence X : = sucn tnat 

Zk = \{i E N such that Xi € X and Xi — a k }\, (3) 

where | • | denotes the cardinality of a set. We write X — 

Mz). 

Each representation (configuration Z or position sequence 
X) has notation advantages and drawbacks. For example, 
given a configuration, transitions like steps, attachments and 
detachments will be described by adding a proper "transition 
sequence". Conversely, to determine the forces the cargo 
particle is subject to the representation it terms of a position 
sequence will be easier to handle. 

Definition 3: We define Z m as the set of all the possible 
configurations with m engaged agents and X m as the set of 
all possible position sequences with m engaged agents. We 
also define 



Z := |J Z v 



(4) 



m— 1 



m— 1 



Lemma 4: The set Z is countable. 

Proof: The proof is left to the reader. ■ 

Lemma 5: For any Z € Z, let X = 4>a{Z) E X be the 
position sequence of the configuration Z with respect to the 



location sequence A. The function <\> A : Z — > X is a one- 
to-one correspondence. Moreover, i^(.Z TO ) = <Y m for any 
to e N. 

Proof: The proof is left to the reader. ■ 

Definition 6: For any a £ Z and Z = {z k } ke z € Z, we 
define the shift operator p Q : Z -> Z 

:= {z fe+Q } fceZ . (5) 
We assume that the motors move towards larger values 
of their positions. With this assumption, we introduce the 
following definition. 

Definition 7: Consider Z = {z k } ke z € Z and let X = 
4>a{Z) — (x\, ...,x m ) be its position sequence with respect 
to the location A. We say that x m is the vanguard agent 
position while X\ is the rearguard agent position of the 
configuration Z (and of the position sequence X). We define 
spread of Z (or spread of X abusing the nomenclature) as 

s{Z) := x m - xi (6) 

and the dimension of Z as 

h(Z) := 1 + sup{fc e Z|z fe 7^ 0} - inf{fc e 7L\z k ^ 0}. 

(7) 

If m = we define s(Z) := and := 1. 

In the model, only certain transitions are allowed, rep- 
resenting a forward step, a motor detachment from the 
microtubule and a motor attachment to the microtubule. 

Definition 8: Define the bi-infinite sequences 

R k s) ■= 4+i - 4 (8) 
R[ d) ■= -5 k (9) 

R k a) := S k (10) 

where S k represents the bi-infinite Kronecker delta sequence 
that is equal to 1 at k and zero anywhere else. 

Definition 9: Consider a time-continuous Markov Process 
defined on the state space Z, with time variable t G R and 
with transition rate from the configuration Z x to Z 2 given 
by X(Z2, Zi), meaning that, given that the system is in the 
configuration Z\ at time, the probability of being in Z 2 at 
time t + At is equal to X(Z 2 , Z x )At (for small At). 
Assume that, for any Z 2 , Z\ such that X{Z 2 , Z\) > 0, there 
exists k € Z such that one of the following conditions is met 

• Z 2 = Z\ + 

• Z 2 = Z\ + R k 

• Z 2 = Z\ + R k . 

We call the transitions sequences R k s \ i?^ and R, a ^ re- 
spectively forward step from location k, detachment from 
location k and attachment to location k. We call such a 
process an Ensemble of Molecular Motors (EMM). 
The cargo transported by the group of molecular motors is 
usually subject to a constant force Fi oa d- In this simplified 
article it is assumed that, after any transition, the cargo 
reaches its equilibrium position instantaneously. We provide 
the following definition. 

Definition 10: Assume that a motor in position x is sub- 
ject to a force equal to \{ x ~ given that the cargo is 



in position x^ £ R. 

Given a position sequence X = {x\, x m } (or equivalently 
its configuration Z = </)T 1 (AT)), an equilibrium position of 
the cargo x^ is any real number satisfying the equilibrium 
condition 

m 

Fi oad = J2x(xi-x {c) )- (ID 
i=i 

Definition 11: We refer to any realization Z(t) of a EMM 
as a configuration path. In an obvious way, any configuration 
path Z(t) can be associate to a position path X(t) = 
<f>A{Z(t)) = (xi(t), £ TO (t)(t)). For any t, we also have 
the number of engaged agents m(t), the vanguard agent with 
position x m (t){t), the rearguard agent x\(t) and, under the 
assumption that it is unique, the cargo equilibrium position 
x^ c \t). 

In [6] a set of specific rules are defined to create a stochas- 
tic simulation algorithm for the motion of molecular motors. 
The rules are the following. The cargo has to attached 
motors that can attach and detach from the microtubule. 
Their linkages behave as elastic springs of elastic constant 
K and rest length l that buckle down with no resistance 
when compressed. It is postulated that motors engaged on the 
microtubule have always a chance of detaching depending 
on the force they are subject to. On the other hand, any 
disengaged motor can attach again, but only at locations 
where its linkage is not subject to a stretch (that is at 
distance not larger than l from the cargo). The probability 
of taking a step depends on the force acting on the specific 
motor and there exists a "stalling force" F s > making 
the step not possible. Finally, the behavior of the system is 
invariant to translations along the microtubule. These rules 
are formalized in the following way 

Definition 12: Given a EMM, we say that it satisfies the 
Kunwar-Veshinin-Xu-Gross (KVXG) rules with parameters 
(m, K,lo,F s ), with to e N and K, Iq, F s real positive, if it 
meets the following properties. Assume a location sequence 
A = { a k}kez\ for any Z\ ^ Z 2 € Z, with position 
sequences X\, X 2 and engaged agents mi, m 2 respectively, 
we have 

• irreversible cargo loss with no engaged motors 

mi = implies X(Z 2 , Z x ) = (12) 

• at most to engaged motors 

m 2 > to implies X(Z 2 , Z\) = (13) 

• spatial invariance (for any aeZ) 

\(Z 2 ,Z 1 ) = X(p a Z 2 ,p a Z 1 ) (14) 

• elastic stretch with constant K and dead zone lo 

( K(l - l ) if I > l 
X(l)~{ if|/|</ (15) 

[ K(l + l ) if I < -lo- 
rn attachment with no stretch 

X(Z 2 , Zx) > and Z 2 -Z x = R k a) 
implies \a k — x^ \ < l (16) 



• stalling force F s 

X (a k - x®) > F s implies \{Z l + R^^Z^ = 0. 

(17) 

Given an EMM with a transition rate function A : Z — >• 
R+ meeting the KVXG rules, define P z (Z,t\Z,t ) as the 
probability of being in the configuration Z — {zk}kez at 
time t > to conditioned to Z(to) = Z. 

Then, given an initial time tg and an initial state Z, for 
t > to, Pz(Z,t\Z,to) satisfies the Master Equation 

^,P z (Z,t\Z,to) = -P z (Z,t\Z,t Q ) ]T KZ',Z) (18) 

Z'ez 

+ HZ,Z')P z (Z',t\Z,t Q ), (19) 

Z'eZ 

that represents the conservation law of the probability mea- 
sure. 

III. Finite dimensional projection and Master 
Equation analysis 

The study of the dynamics of the motors (represented 
either in the configuration space or in the position sequence 
space) is not easy to handle since the entities involved are 
not stationary and the state-space is infinite-dimensional. The 
following theoretical result is helpful for transforming the 
problem into a finite dimensional one. 

Theorem 13: Consider an EMM satisfying the KVXG 
rules with parameters (m, K, Iq, F s ). For any configuration 
path Z[t) let s(Z(t)) be its associated spread for any t. 
Define 



s (max) ._ max 



mF, - Fi 



load 



F, 



K 



load 

~K~ 



2l . (20) 



For any S > s( moa: ) and for any t > t , we have that if 

s(Z) < S 



Pr{s(Z(t)) < S | Z(to) = Z} = 1. 



(21) 



where Pr{E\C} is the conditional probability of an event 
E given the condition C. 

Proof: See the appendix ■ 



Theorem 1 3 enables us to study many properties of an EMM 
using a formulation that has finite dimension and deals with 
stationary quantities. To this aim, we introduce the following 
definitions. 

Definition 14: Given an EMM satisfying the KVXG rules 
with parameters (m, K, l , F s ) subject to load Fi oac [ and a 
location sequence A with step size d s , we define 

" „(mas) " 

(22) 



as the regular dimension of the ensemble. 

Given a configuration Z = {zk}kez, we also define its 

ensemble representation as the n-vector 



Q — ( z fci i ■■•^fci+n-l)"' 



(23) 



where k\ := inf{fc e Z\z k ^ 0} and we write Q := U^(Z). 
If fei = -oo, Q is (0,0,...,0) T . 



The ensemble representation of Z is the projection of the 
configuration on a suitable n-dimensional space. Theorem 
13 guarantees that the ensemble representation Q(t) of the 
configuration Z(t) contains all the non-zero entries of Z(t) 
with probability 1 for all t > to if the dimension of Z(to) 
is not greater than n. In this scenario, while Q(t) does not 
contain the information about the absolute positions of the 
motors, it still contains all the information about all their 
relative positions, so it can be used as a tool to study the 
dynamics of their mutual interactions. 
Assume that the starting configuration Z(to) = Z has 
dimension not exceeding the regular one n. For t > to the 
spread is not exceeding n and the transition probabilities 
of Z{t) are known. Despite the projection, it is possible 
to derive the transition probabilities of Q(t). This allows 
to define the dynamics of a Markov system with a finite 
dimension. The following properties hold. 

Lemma 15: Consider an EMM with regular dimension n. 
Let s(Z) < n and P Q {Q,t) be the probability that Ii {e \Z) 
is equal to Q at time t. Then, for Q ^ and for any t > t 

P Q (Q,t\Z,t ) = Y, P z{p aZ ^M (24) 

Proof: The proof is left to the reader. ■ 
Theorem 16: Consider an EMM with regular dimension 
n and let s(t) be the spread at time t. Define 



Aq(Q',Q):= J2 KZ',Z). 



(25) 



n( e )(z') = Q' 
s(Z)<n 



For any t > to and for any Z such that s(Z) < n, it holds 
that 



d_ 

dt 



-P Q (Q,t\Z,t ) = -P Q (Q,t\Z,t ) HQ', Q) (26) 

Q'eQ 

+ Yl HQ,Q')P Q (Q',t\z,t ) (27) 

Q'eQ 

Proof: In order to simplify the notation we omit the 
initial condition: P Q {Q,t) = P Q (Q,t\Z,to) and P z (Z,t) = 
P z (Z,t{Z,t ). We have that Q = if and only if Z = 0. 
In such a case, X(Z' ', Z) = for any Z' and X(Z, Z') = 
for any Z' such that U^(Z') ^ Q 1 := (1, 0, 0). Observe 
that 

• Aq(Q',0) = for any Q' 

• A Q (0,Q') = 0forany QVQi 
. P Z (Z',t) = if s(Z') > n. 

Then, for Q = we obtain 



^(0,0 = ^(0,*) = 

= mz')Pz(Z , ,t) = \ Q (0,Q 1 )P Q (0,Q 1 ). 

n< c >(z')=Qi 
s(Z)<n 



9 Pa(Q,t) = §-P Q (n {e) (z),t) = 53 l f Pz(p a z,t) = 



For Q 7^ 0, we have 

d d 

a6Z 

-EE Kz',P a z)Pz{ P a z 1 t) 

+ E E Kp a Z,Z')P z {Z',t) = 

aez Z'e2 

= -E E A(p^',p tt Z)F z (p«Z,i) 

+ E E X(p a Z,p a Z')P z {^Z',t) 
aez Z'e2 

=-53 a^zjx;^^*) 

+ ^ A(Z,Z')^Pz(p^',t) = 

= -Pg(Q,t) ^ A(Z',Z) 
Z'ez 

+ ^ A(z,z')P Q (n( e )(z')) = 

s(Z')<n 

= -P Q (Q,t) 53 A Q (Q',Q) 
Q'eQ 

+ E aq(q,q')W) 

Q'eQ 

■ 

Lemma 17: Given an EMM satisfying the KVXG rules 
with parameters (m, K, lo,F s ) with regular dimension n, the 
number of possible ensemble representations is given by 



N= 1 + 



\ " 



(n + m- 2)! 



^ fn-lWrn- 1)!' 

m=l v ; v ' 

Proof: The proof is left to the reader. 



(28) 



Enumerate all the possible ensemble representations 



Qi,-~,Qn with N given by (28 1 and define Pj(t) as the 



probability of having the ensemble representation Qj at time 
i given the initial condition Z such that s(Z) < n. By 



projecting Equation (18i onto the set Q := {Qi, Qn} 
we obtain 

d 
~dt 



;P Q (Q,t) = - 53 X Q (Q',Q)P Q (Q,t) (29) 



Q'eQ(+)(Q) 



Q'6QC-)(Q) 



*Q(Q,Q')PQ(Q',t), (30) 



The dynamics of the probabilities on the space of the 
ensemble representations is given by an equation of the form 



— P = AP 
dt 



(31) 



where P := (Pi, Pn) t , the off-diagonal element of A, 
dji is the rate of the transition from Qi to Qj and the sum 
of all the rows of A is zero. In this way the dynamics of 
the infinite dimensional Markov system has been projected 
onto the dynamics of a finite dimensional one. In the section 
we show how this projection provides no loss of information 



when used to obtain certain quantities of interest such as the 
expected velocity of the cargo and the expected runlength. 

IV. Determining the transition probabilities 

FROM SINGLE MOTOR OBSERVATIONS 

In this section we determine the transition rates X(Z' , Z) 
for the infinite dimensional stochastic model. The transition 
rates \q(Q',Q) necessary to define the finite dimensional 
one can be computed from the transition rates X(Z',Z). 
Following [6], it is possible to derive the probability rates 
of simple events X(Z', Z) by considering the behavior of a 
single motor. 

Rate for stepping transitions 

Every motor moves converting ATP into kinetic energy 
according to the chemical equation 

K+ATP t± K ATP ^ ' K+ADP+Pi+energy. (32) 

k °tf 

Following [9], a Michaelis-Menten dynamics leads to a 
hydrolisis rate equal to k cat [ATP]/([ATP]+k m ). It is also 
assumed that the free head of the motor successfully binds to 
the microtubule location with probability rate (or efficiency) 
E. Then, the transition rate X step of stepping and the average 
velocity V of a single motor are given by 



_ k cat [ATP] 
step ~ [ATP] + k m £ 



V 



Pstepd 



k cat [ATP] 

[ATP] + krr. 



de 



(33) 
(34) 



where d is the fixed size of the protein step. The force F that 
the cargo exerts on the motor is assumed positive when it 
opposes the motor motion. It is assumed that when it exceeds 
a certain force F$ it causes the motor to stall. Following [6], 
it is assumed to affect the motor dynamics by changing the 
probability e of binding with the microtubule, following the 
relation 



1 if F < 

--(F) = { 1 (/) 2 if 0<F<F o 
otherwise 



(35) 



where F s is the stalling force of the KVXG rules. In Figure [2] 
it is represented how the efficiency e(F) changes in different 
situations: when the motor is lagging behind the cargo, when 
it is at a distance where the linkage is not stretched and when 
it is actively pulling the cargo. 

The linkage is modeled as an elastic spring of elastic 
constant K and rest length l that buckles down with no 
resistance when compressed. The relation between the force 
applied to the motor and the linkage length / is given by 



K(l + l Q ) if I < -l a 
F(l) = { if|Z|<Z 
K(l - l ) if I > l . 



(36) 



In [6] it is assumed that the force F also influences the 




Fig. 2. Graphical representation of the "efficiency" e(F) and its depen- 
dence on the force F applied to the motor 



K eI (l - l a ) 



la I 



kinetics of the ATP hydrolysis. In particular it is assumed 
that k ff increases with increasing F 



k tt - kn ttP Fd '/ K l> T 



(37) 



where fco / / is the backward reaction rate of the hydrolysis 
when F = 0, Kb is the Boltzmann constant, T is the 
temperature and di is a parameter that can be experimentally 
determined. Thus, the transition rate for a step at location a/, 
is given by 



X(Z + R{ s) ,Z) = z k 



k ca t[ATP] 



ATP] 



kcat[ATP] 



V Fa 



k on +k off (F) 

kaat 

2 



[ATP] 



kaat 



e{F) (38) 



(39) 



Rate for detachment transitions: From [Schnitzer et al. 
2000], the processivity L is 



-F8,K h T 



L = 



d[ATP]Ae 
[ATP] + B{1 + A)e- F ^ K » T '' 



(40) 



where A, B and 6i are parameters that can be experimentally 
determined. Since the processivity represents how far a motor 
can move before detaching from the microtubule on average, 
we find a relation between the probability of stepping and 
the probability of detachment. 



step 



L 



[ATP)Ae- FSlKbT 



Pdetach(F) d [ATP] + B(l + A)e-^ b T ■ 
Thus, so long as F < Fq, 



(41) 



[ATP]+B{l + A)e-™*> T 

^detachVf ) - \ATP]Ae- F5lKbT r stepV-T )■ 



(42) 



When F > F s , in [6] a constant detachment rate is assumed 
Pdetach(F) = Pback = 2s" 1 . Thus, we have the following 



transition rates 



\(Z + R[ d \z) 



[ATP]+B(l+A)e- Xia i'-^ '>*> 
[ATP]Ae-x( a k~^ c) hK b i 

if xK - x(c) < F s 

^k-Pback 

if xK - ^ (c) > F s 



'-\{Z + R { ° 



(43) 



Probability of attachment: Experimentally, it is found in 
[6] that the probability of a motor to attach to the microtubule 
is P at t — 5s _1 . If the motor is linked to the cargo, it is 
assumed that is attaches to the microtubule without stretching 
its linkage. Thus, the only admissible locations of attachment 
are the locations at a distance from the cargo that is less than 
Iq. They are also assumed all equally likely. 

V. Probabilistic Analysis 

We have shown that the length of all the possible strings 
can be uniformly bounded. As a consequence, apart from 
possible spatial shifts along the microtubule, the number of 
possible arrangements of the motors is finite. 

Using the finite dimensional Markov Model derived in 
the previous sections, it is possible to derive meaningful 
characteristics of the system in an exact manner. 

A. Transient analysis 

In Figure [3] we present the dynamics of P(t) in the case 
m = 3 for two different values of Fi oa d- We use 4 time 
snapshots to describe how P(t) changes with time. In the first 
row of Figure [3] we present the probability density function 
of the configurations when the cargo is subject to a low 
load (Fi oa d = 0.0002niV). The x and the y axis represent 
the distance of the two leading motors from the rearguard 
motor. The distance is expressed in number of locations 
on the microtubule. The z axis represents the probability 
of having a motor at distance x and a motor a distance y 
from the rearguard motor. The probabilities of configurations 
with less than three motors engaged are not reported be- 
cause that would make the visualization unnecessarily more 
complicated. The four snapshots are taken at time Osec, 
0.2sec, 0.4sec, 0.6sec from left to right. In the second row 
of Figure [3] we report the results in the case of a high 
load {Fi oa d = 0-008nN) with an identical initial probability 
density function. Notice how the change in the parameter 
Fioad leads to two different dynamics of the probability 
density function: under a low load the motors tend to spread 
out while under a high load certain configurations tend to be 
more likely. 

1) Average runlength: The detailed knowledge of the 
probability vector P(t) can be used to derive information 
about numerous aspects of the behavior of the motor ensem- 
ble and study the variation of certain quantities as a function 
of different parameters. An important quantity that can also 
be experimentally measured in experiments is the expected 
runlength of the motors, that is the average length travelled 
by the cargo before being lost. The average runlength is 






immediately computed from P(t) 



Average Runlength = E 



+ °° dx^ 
dt 



dt 



(44) 



where E[-] is the mean operator and x^ is the cargo 
position. This same quantity as a function of the load has 
already been computed in [6] using Monte Carlo methods in 
two scenarios. In one scenario (Model A) the authors neglect 
the effect of thermal noise and in another scenario (Model 
B) they provide a dynamic model for the brownian motion 
of the cargo. 

We report our results in Figure |4] for two reasons. The first 
reason is for comparison. By forcing the variance of the 
cargo position to zero we obtain a model that is equivalent 
to "Model A" in [6]. The only difference is that in our case 
we are computing the average runlength from the probability 
density function without using a simulation approach. Thus, 
forcing the variance of the cargo position to zero, the results 
must match. We report the results of our computation in 
Figure [4] a using both a coarse grid (solid lines) and a fine 
grid (dashed lines) for the load force. If we compare the 
results we computed with the results in [6] at the same 
points (that is using the coarse grid), the runlength-load 
curves overlap almost perfectly. The second reason is to 
highlight a possible weakness of the model adopted for the 
detachment rate. In the finer scale, we can notice peaks in 
the runlength curves that were not evident before. These 
peaks appear in correspondence of loads that are multiple 



of the stalling force Fq and can be explained as "artifacts" 
of the model. They can be explained in the following way. 
Let us consider the case m = 2 for simplicity. When only 
one motor is engaged (m(t) — 1) and Fi oai i is close to (but 
less than) the stalling force Fq, the probability of detachment 
becomes small, and the loss of the cargo becomes unlikely. 
In the meantime, the disengaged motor has the opportunity of 
attaching to microtubule, catching up with the leading motor 
and moving the cargo a little further. For m > 2, similar 
arguments can be used. This mechanisms explains how, in 
the mathematical model, by neglecting the brownian motion 
of the cargo, the expected runlength tends to increase while 
the load approaches multiple values of the stalling force. 
When the brownian motion of the cargo is taken into account 
(see Figure [4] b) those peaks are smoothed down, but do 
not disappear completely. Thus, the pronounced peaks in 
Figure [4] a and Figure [4] b are due to the inaccuracy of the 
equation ( |38] > for F ~ F s . Still, the insights given by the 
model could have a physical meaning and could be used to 
explain possible mild breaks of monotonicity in the runlength 
curves for multiple motors. 

B. Steady state characteristics 

1) Average number of engaged motors: In Figure [5] a, 
the probability distribution P(t) has been used to determine 
the probability of having a given number of motors engaged 
on the microtubule at time t. We notice that the probability 
of cargo loss increases with time (and it can be proved 



(a) (b) 

Fig. 4. Average runlength neglecting thermal noise on the cargo (a) and considering it (b). 





(a) (b) 

Fig. 5. Histogram of the number of motors (left) and histogram of the number of motors assuming that the cargo has not been lost, that is m(t) > 0. 





(b) 



Fig. 6. 



converges to 1 for t — » +oo). Thus, it is not possible to 
define a steady state for the system: the eventual cargo loss 
is the actual steady state. However, in Figure [5] b we report 
the conditional probability density function of the number 
of motors engaged on the microtubule at time t, given that 
at least one motor is engaged. We observe that using this 
conditional probability density function a steady state can be 
defined. In the case of m = 3, the steady state conditional 
probability is reported in the histograms of Figure [6] a and in 
Figure [6] b for Fload = 0.0002niV and Fload = 0.008niV 
respectively. As in Figure [3] the rearguard motor is taken as a 
reference point, the x and y axis represent the distance of the 
other two motors from the rearguard motor and the z axis 
represents the probability. It can be observed that under a 
low load the motors tend to be spread and distant from each 
other, while under a high load they tend to be clusterized. 
Indeed, the pronounced probability peak around the point 
x = 0, y = indicates that the three motors are not far 
from each other. It is interesting to notice that there is a 
significantly higher probability along the diagonal line x = y. 
This can be interpreted as follows. The three motors tend to 



be close together, but when they are not it is more likely 
to find the two leading motors together while the rearguard 
motor is lagging behind. 

2) Average velocity: The definition of the steady state 
conditional probability can be used to determine the average 
steady state velocity of the motor ensemble. The results are 
reported in Figure [JJ a for the case where thermal fluctuations 
are neglected and in Figure [7J b for the more accurate 
model. Results in Figure [7J a are again for comparison 
and verification since, in the noiseless scenario, our model 
is equivalent to Model A in [6]. Results qualitatively and 
quantitatively non-dissimilar are reported in Figure [JJ b for 
the noisy scenario. 

C. Detection of rare events 

The possibility of determining exactly the probability 
distribution of the different motor configurations allows for 
the detection of rare events. It is possible, for example, to 
determine the probability of the different steps sizes for an 
ensemble of 2 motors as represented in Figure [8] As it can 
be noted, there are two prominent peaks corresponding to 
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Fig. 




(b) 



Fig. 7. 



8 nm and 4 nm. These peaks corresponds respectively to 
the case where there is one active motor before and after 
the step and to the case where there are two active motors 
before and after the step. There are also different predicted 
step sizes at about 2 nm, 6 nm and around 4 nm. They 
correspond to situations where there is a different number of 
active motors before and after the step and they can be mostly 
considered as artifacts of the model. However, we also find 
a small probability, but different from zero, of steps larger 
than 8 nm at about 11 nm. This does not seem possible, 
since the each motor can advance only 8 nm at the time. 
Indeed, this can occur when the rearguard motor is pulling 
back the cargo and detaches from the microtubule. While 
this situation is not very likely in a bead-assay experiment 
that our simulations are trying to portray, it is quite likely 
that two motors could pull in two opposite directions in a 
gliding assay experiment. Thus, steps larger that 8 nm as 
those described in [10] could be originated by a mechanism 
of this kind. The probability distribution of the step size for 
an ensemble of 3 motors is reported in Figure [8] where the 
steps longer than 8 nm have an identical interpretation. 

Appendix 

A. Preliminary Lemmas 

Lemma 18: If Fi oa d > and m > 0, then the equilib- 
rium solution x^ to (111 for any position sequence X = 
(xi, x m ) and for \ as given by ( 15 1 exists, is unique and 
x rn ^* ^ ^- 

Proof: The proof is left to the reader. ■ 

Lemma 19: Consider an EMM meeting the KVXG rules 
with parameters (m, K,Iq,F s ). Assume that Fi oa d > 0, and 
that, for a configuration Z, the number of engaged agents is 
m > such that the cargo position x^ is uniquely defined. 
Also assume X = 4>a{Z) = (xi,...,x m ) for a location 



sequence A and that 

s(Z) := x m 



xi > h 2l . (45) 

A 

Then the vanguard motor is stalled, that is F s < x( x m ~ 
x^ c '). If the vanguard motor is in the location <jfc m it imme- 
diately follows that, for k > k m , \{Z + R^(k\Z) = 0. 

Proof: Without any loss of generality, assume x^ = 0, 
since the reference point on the microtubule is arbitrary. By 



contradiction, and by using Lemma (18 1, assume K{x r 
Iq) < F s . If x\ < —Iq, then we have 



F, 



load 



K(x x + Jo) = F t 



oad 



< 



< (rn - l)K(x m - Iq) < -K(x m - l ) + mF s 
This implies 

mF s - Fi , 



Xl < 



K 



2h 



(46) 
(47) 

(48) 



that is a contradiction. In the case x\ > —Iq, we have 

m 

Fioad = < mF s < mF s + K(x m - l ). (49) 

i=\ 

implying 

, , mF s Fi oa d 07 /cn\ 

x m ~ xi < x m + l < — 1- 2Z , (50) 

A 

that is again a contradiction. ■ 
Lemma 20: Consider an EMM meeting the KVXG rules 
with parameters (m, K, Iq , F s ) . Assume that Fi oa d > 0, and 
that, for a configuration Z, the number of engaged agents is 
to > such that the cargo position x^ is uniquely defined. 
Assume that s(Z) — x m — X\ < q, with 

Fioad 



q> 



K 



2l , 



(51) 



and that the next transition occurring is an attachment at 
position k, that is Z' = Z + #0) (k). Then, if X(Z', Z) > 

s(Z') < q. (52) 
Proof: Without any loss of generality, assume x^ = 0, 
since the reference point on the microtubule is arbitrary. Let 
X' = x' m+1 ) — (Pa(Z') If xi < —lo, the transition 

does not affect the value of the distance because the newly 
attached motor will be in the range [—Iq, Iq] in order to meet 
the KVXG rules. Indeed, we have 

s(Z') = x' m+l - x[ = x rn -xx = s(Z) < q (53) 

If x\ > — Iq, then x' x > — Iq- Moreover, 



Fload = x( x i) > K{x m -l ), 



i=l 



(54) 



implies 



x' m+l - x[ < x m + l < + 2l < q. (55) 



B. Proof of Theorem \13\ 

Proof: Define, for t > to, 

< <f>(t) := Pr{s{t) > S\Z(t ) = Z}= (56) 
= £ P z (Z,t\Z,t ). (57) 

s(Z)>S 

At time to we have that <f>(to) — 0. By contradiction assume 
that there exists a time t\ > to such that 4>{t\) > 0. We must 
also have that, for some t E (to,ti), 



^(t) > o. 



(58) 



From the master equation 





dt 



E I" E HZ',Z)P z (Z,t\Z,t ) 
s(z)>s I z'ez 



X(Z,Z')P z (Z',t\Z,t ) 



Z'ez ) 
= E E {-KZ',Z)P z (Z,t\Z,t ) + 

s(Z)>S s(Z')<S 

+\(Z,Z')P z (Z',t\Z,t )} < 
< E E {HZ,Z')P z (Z\t\Z,t )} 

s(Z)>S s(Z')<S 

Now we prove that, if s(Z') < S and s{Z) > S, we 
necessarily have \(Z, Z') — 0. 

The only possible transitions are detachment, attachment and 
step. 

The detachment transition can only reduce the configuration 
distance, so Z - Z' ^ (k) for any k. 



If the transition is an attachment, by Lemma 20 we obtain 
the assertion. 

If the transition is a step the configuration distance increases 
only if the vanguard motor steps, and the increase is by d s . 



If s(Z') < S — d s , then s(Z) < S and this is not a case 
we are considering because we have assumed s(Z) > S. 
Conversely, if s(Z') > S — d s , by Lemma 19 we have that 



the vanguard motor can not step, so such a transition is not 
possible. Then we obtain as a contradiction that, for t > to, 

J^(t) < 0. (59) 
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